home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libconv / TRY / ornl_diir1d.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  4.7 KB  |  137 lines
  1.       subroutine ornl_diir1d(    in_put, incinp, i0_inp, n_inp,
  2.      $                iirfil, inciir, i0_iir, n_iir,
  3.      $                output, incout, i0_out, n_out)
  4.       double precision in_put(*), iirfil(*), output(*)
  5.       integer incinp, i0_inp, n_inp
  6.       integer inciir, i0_iir, n_iir
  7.       integer incout, i0_out, n_out
  8. c-----------------------------------------------------------------------------
  9. c
  10. c   dconv1d  performs a 1D recursive convolution in the time domain :
  11. c    O(j) = (1/F(0)) * { I(j) - Sum[ O(i) * F(j-i) ] }, for i=1,...,(j-1)
  12. c
  13. c-----------------------------------------------------------------------------
  14. c
  15. c   PARAMETERS:
  16. c
  17. c    in_put:    Pointer to FIRST sample of sequence "in_put"
  18. c    incinp:    Increment between two successive values of "in_put"
  19. c    i0_inp: Index of the first element of "in_put" 
  20. c    n_inp:    Number of samples of "in_put"
  21. c
  22. c    iirfil:    Pointer to FIRST sample of sequence "iirfil"
  23. c    inciir:    Increment between two successive values of "iirfil"
  24. c    i0_iir: Index of the first element of "iirfil" 
  25. c    i0_iir:    Number of samples of "iirfil"
  26. c
  27. c    output:    Pointer to FIRST sample of sequence "output"
  28. c    incout:    Increment between two successive values of "output"
  29. c    i0_out: Index of the first element of "output" 
  30. c    n_out:    Number of samples of "output"
  31. c
  32. c-----------------------------------------------------------------------------
  33. c
  34. c   PLEASE NOTE:
  35. c    1-    Please note that the array pointers must all point to the first 
  36. c    element of the array "i0_inp", "i0_iir" and "i0_out".
  37. c     If "in_put" for example is defined as
  38. c        dimension in_put(-25:45)
  39. c    Then "dconv1d" must be called with the following parameters
  40. c        call dconv1d( in_put(-25),1,-25,45, ... )
  41. c
  42. c    2-    There isn't any checking done to see if the IIR filter is minimum phase
  43. c    or not. In the case of a non-minimum-phase filter, the computation is
  44. c    numerically unstable, OVERFLOW may results.
  45. c
  46. c-----------------------------------------------------------------------------
  47. c
  48. c   USAGE:
  49. c    This module compute the result of the convolution in the "output" range
  50. c    padding with zeroes when needed.
  51. c
  52. c    With some precautions it is possible that 
  53. c        the Output sequence   OVERWRITE   the Input one.
  54. c    Then (if i0_out == i0_inp), it is necessary that   i0_iir =< 0
  55. c    Said in DSP jargon: "iirfil" must be ANTI-CAUSAL
  56. c
  57. c    In theory, an input sequence of "n_inp" samples starting at time "i0_inp", 
  58. c    filtered by a sequence of "n_fir" samples starting at time "i0_fir",
  59. c    will result in a new signal of infinite length starting at time (i0_inp + i0_fir). 
  60. c    He we just compute here the values that fall in that range and zero the rest.
  61. c
  62. c     For example when filtering a sequence of N samples, with a
  63. c    filter of m samples. If one wants only to compute the first N resulting
  64. c    samples, the following call can be used:
  65. c        call diir1d( f, 0, 1, N,  g, 0, 1, m,  h, 0, 1, N)
  66. c
  67. c-----------------------------------------------------------------------------
  68. c
  69. c   This Fortran Subroutine written by
  70. c    Jean-Pierre Panziera
  71. c    Silicon Graphics Inc
  72. c    September 20, 1991
  73. c
  74. c-----------------------------------------------------------------------------
  75.  
  76.       double precision tmp, zero, one, alpha
  77.       parameter (zero = 0.0, one = 1.0)
  78.       integer i1_inp, i1_iir, i1_out, k0, k1
  79.       integer i, j, kout, kinp, kiir, n0, n1, jmax
  80. c-----------------------------------------------------------------------------
  81. c
  82. c   First Check Parameters
  83. c
  84. c-----------------------------------------------------------------------------
  85.       if( iirfil(1) .eq. zero ) then
  86.     print *, 'In - diir1d - Error - First sample of IIR filter must be != zero'
  87.     stop
  88.       end if
  89.  
  90.     i1_inp = i0_inp + n_inp - 1
  91.     i1_iir = i0_iir + n_iir - 1
  92.     i1_out = i0_out + n_out - 1
  93.     k0 = i0_inp - i0_iir
  94.  
  95.     n0 = i0_out
  96.     if( n0 .lt. k0)        n0 = k0
  97.  
  98.     alpha = one / iirfil(1)
  99. c-----------------------------------------------------------------------------
  100. c
  101. c   Zero the output at the begining of the computated window
  102. c
  103. c-----------------------------------------------------------------------------
  104.     kout = 1
  105.     jmax = min ( i1_out, n0-1)
  106.     do j = i0_out, jmax
  107.         output(kout) = zero
  108.         kout = kout + incout
  109.     end do
  110. c-----------------------------------------------------------------------------
  111. c
  112. c    Compute the convolution
  113. c
  114. c-----------------------------------------------------------------------------
  115.     do j = j, i1_out
  116.         if( j .le. (i1_inp-i0_iir) ) then
  117.         kinp = 1 + (j+i0_iir-i0_inp) * incinp
  118.         tmp = in_put(kinp)
  119.         else
  120.         tmp = zero
  121.         endif
  122.         k0 = max( j-(n_iir-1), n0)
  123.         k1 = j-1
  124.         kout = 1+ (k0-i0_out) * incout
  125.         kiir = 1+ (j-k0)*inciir
  126.         do i = k0, k1
  127.         tmp = tmp - output(kout) * iirfil(kiir)
  128.         kout = kout + incout
  129.         kiir = kiir - inciir
  130.         end do
  131.         kout = 1 + (j-i0_out) * incout
  132.         output(kout) = alpha * tmp
  133.     end do
  134.  
  135.       return
  136.       end
  137.